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I. INTRODUCTION 

Liquid metals are complex binary fluids consisting of ions in a sea of conduction electrons. While the ions can usually 
be treated classically, the electrons are typically degenerate and must be treated quantum-mechanically. Liquids are 
differentiated from gases by non-trivial structure at the level of two-body correlation functions; they are generally 
close in density to solid phases. For two-component systems these correlation functions are defined in /c-space as: 

< p a (k)p>(-k) > 1/2 
OupW = (jv Np) 1 / 2 ( 7V «- /N W °k,o- (l-l) 

The S a f3(k) are referred to as static structure factors and the operator 

N a 

p Q (k)=$> 4kriQ ' ( L2 ) 



is the Fourier transform of the one-particle density operator of component a. The indices a and (3 refer to ions, (/), 
or valence electrons, (e). The structure factors S a p(k) can be related to the so-called radial distribution functions 
g a p{r) by: 



S a p(k) = 5 a p + { Pa ppY 2 / dre lk r [g a p(r) - 1] , (1.3) 
Jv 

where the pi are the homogeneous average densities. 

The determination of the ion-ion structure factor Sn{k) and the electron-electron structure factor S ee (k) are 
interesting problems in their own right (one largely quantum mechanical, the other largely classical), and have been 
the focus of much research: the Sjj(k) because of their experimental accessibility; the S ee (k) (with the ions usually 
smeared into a rigid neutralising background) because of the importance of the electron fluidtl. In contrast, the 
electron-ion structure factor S e i(k) has received considerably less attention, partially because it is hard to measure, 
partially because its exact physical relevance remains largely unexplored and unknown, and partially because it 
includes both the physics of the ions and the physics of the electrons, each of which is traditionally treated with its 
own set of theoretical techniques. 

One of the simplest-ways to treat the valence electrons in a liquid metal is in a linear response formalism using a 
local pseudo-potentialo. In fact, linear response has been shown to be much more accurate than one would naively 
expect, a result which stems in part from a recently discovered interference effect between an atomic lengthscale, 
the inverse ionic length, and an electronic lengthscale, twice the Fermi wave vector 2fcpLl. This interference effect 
significantly reduces the magnitude of the non-linear response terms at the normal densities of most liquid metals. 
The electron-ion correlations emerge when the induced linear response electron density is combined with standard 
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liquid state techniques to treat the ionsmJ. This approach is easy to implement, can in some cases be remarkably 
accurate, and can explain the qualitative trends in the shape of the electron- ion structure factor S e i(k) for metallic 
liquids across the periodic tablet! The main obstacles to higher accuracy lie in the uncertainty over the exact (local) 
pseudo-potential, especially when non-local effects are important!], and f\f^i in the neglect of non-linear electron 
response and of ion-ion correlation effects on the induced electron densities&t£l. ._. 

The development of ab initio simulation techniques based on Density Functional Theory (ftFT) for the electronsliil, 
and molecular-dynamics on the adiabatic electronic potential energy surface for the ionsL3, provide probably the 
most accurate and well-tested approach to electron-ion structure. However, the drawback of these methods is their 
computational cost; in practice only relatively smalLsystem sizes can be investigated and so far only results for Mg 
and Bi electeop-ion correlations have been publishecO. The related Orbital- Free ab initio molecular dynamics method 
(OF-AIMD)tLa allows larger system sizes and significantly longer simulation times, and has been successfully applied 
to the electron-ion correlations of Li, Na, Mg, and Alli3'La, but the computational cost is still rather large. 

An alternative approach is the Quantum Hypernetted Chain (QHNC) method of Chiharalla, which self-consistently 
combines integral equation techniques from the theory of simple liquids with a Kohn-Sham type treatment for the 
electrons. The QHNC treats the electrons and ions on essentially equal footing, does not require a pseudo-potential 
approximation, and is computationally relatively cheap. Ion-ion and electron-ion correlations emerge in the thermo- 
dynamic limit - there are no finite size effects. 

In section II, we derive the basic form of the QHNC approximation by focusing first on the exact Quantum 
Ornstcin Zcrnikc (QOZ) equations in section II. A, and then outlining the approximations needed to derive the QHNC 
approximation in section II. B. The numerical implementation of the QHNC is detailed in the appendix. 

In section III, we describe the ion-ion and electron-ion correlations that emerge from the QHNC for our set of 
metals: Li, Be, Na, Mg, Al, K, Ca, and Ga. 

Even though the valence electron distributions are changed in a bonded environment, X-ray scattering off liquid 
metals has traditionally been interpreted with a free-atom form factor. In section IV, we describe the difference 
between extracting ion-ion structure in X-ray scattering with a free-atom form factor and extracting ion-ion structure 
with a metallic-atom form factor. The effects of bonding on the coherent X-ray scattering intensity may be measured by 
comparing X-ray and neutron scattering determinations of the ion- ion structure factor Sn(k). However, experiments 
and theory have yet to converge on this issue. 

Finally, we present some concluding remarks in section V, and describe some details related to the numerical 
implementation of the QHNC in the Appendix. 



II. QUANTUM HYPERNETTED CHAIN APPROXIMATION (QHNC) 



A. Quantum Ornstein Zernike Relations 



The Quantum Ornstein Zernike Relatiq«ft-|(QOZ) for a two-component system are most naturally derived in the 
context of density functional theory (DFT)Bt3. First we define thaJIclmholtz free-energy for a two-component system, 
which is a unique functional of the two one-body density profileala: 

F[p x ,p 2 ] = F?\p 1 ]+Ff\p 2 ]+F a '\p 1 ,p 2 ]. (2.1) 

The functional is split in the usual way between ideal (non-interacting) and excess (interacting) parts. We then 
introduce the external potential field: 

*a(r) =/ia-0a(r), (2.2) 

which is defined in terms of the chemical potential fi a of species a and the external potential, <f> a (r) which acts on 
species a only. A Legendre transform with respect to these external fields obtains the grand potential: 

0[*l,^ 2 ] = F[ Pl ,p 2 ] + J drpi(r)*i(r) + J drp 2 (r)V 2 (r), (2.3) 

which is in turn a unique functional of the two external potential fields and 

The first two functional derivatives of the Helmholtz free-energy functional w.r.t. the one-particle densities are: 

T~~"~7~T = *«(r), (2.4) 
5p a (r) 



2 



and 



S 2 F 



=X Q;3 (r,r ) 



8p a (r)6p/3(r') 8pp(r' 

The first two derivatives of the grand potential functional w.r.t. the external potential field are: 

59, 



and 



s 2 n S Pa (r 



Xa^(r,r'), 



(2.5) 



(2.6) 



(2.7) 



which defines the susceptibility matrix or matrix of the linear response functions x«/3(r, r'). Thus the two 2nd functional 
derivatives are each cithers' functional inverse, a natural consequence of having two generating functionals linked by 
a Legendre transforrrO. 

The direct correlation functions C a p{r, r') aLan arbitrary two-component mixture are defined in the usual way as 
functional derivatives of the excess free energyc3: 



-jC a0 (r,r>) 



S 2 F e 



6p a (r)Spp(r')' 



(2.8) 



If we then define (xSl) 1 as ^ ne mverse susceptibility matrix of the ideal system, we arrive, by combining equations 



( [2.l| ), (2.5) and (|2.S|), at the following relationship between two 2-by-2 matrices: 

-1 



■C a f3 = (Xafi) 1 - {X%)~ 



(2.9) 



the Quantum Ornstein Zernike Relations (QOZ). They follow from simple properties of the two free-energy 
functionals and in this form they are valid for any two-component inhomogeneous quantum system (the generalis ation 
to more than two components is straightforward). In the homogeneous limit the direct correlation functions of Eq. (2.S) 
reduce to the usual direct correlation functions first introduced by Ornstein and Ze 
that we will be using them throughout the rest of this paper. 

For classical species, the fluctuation-dissipation theorem relates the response functions to density-density correlation 
functions^: 



and it is in this sense 



lim Xa/3(k, 0) = -f3(p a pp) 1/2 S af 3(k), 



(2.10) 



written here for a homogeneous system and in terms of the structure factors defined in Eq. (1.1). For a liquid metal, 
where t he i ons are viewed as classical but the electrons quantum-mecha nical , inverting the matrix in the QOZ relations 
of Eq. (2.9), and applying the fluctuation-dissipation theorem of Eq. ( 2.10 ) for xn(k) and Xei(k) results in: 

S n (k) - [1 + xlPJ(k)C ee (k)/(3}/D(k) 

Sel(k) = 



P -X^{k){C eI (k)/(3)/D(k) 

Pe 



Xee{k)=X^{k)(l-piCu{k))/D{k) 

D(k) = [l- Pl C n {k)\ [l + X { °\k)C ee (k)/p] + PlX { ° ) (k)\C eI (k)\ 2 /(3, (2.11) 

where Xee'(fc) is the well known Lindhard functio?^, the response function of the non-interacting electron gas. In the 
limit thsdrJioth species are classical, the QOZ relations reduce to the usual classical two-component Ornstein-Zernike 
relationsE3. ^ 

The QOZ relations for a liquid metal appear to have been first derived by ChiharaEl Later Ichimaru et. alx3 
derived similar equations from a two-component linear res pon se formulation. The two formulations are equivalent 
if the definitions of the direct correlation functions of Eq. ( |2.8| ) are linked in the usual way to the local field factors 
G a0 (k): 

C a0 (k) 



-V af3 (k){l - G aP (k)), 



(2.12) 



where V a p(k) is the direct potential between species. 
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B. Quantum Hypernetted Chain Approximation 



To solve the QOZ relations for aJiquid metal we recast them into a slightly different form using two stepstl: The 
first step is to use the Percus trickEJ to relate the homogeneous two-body pair-correlation functions to the one-body 
inhomogeneous density around one particle fixed at the origin. For the electron-ion pair-correlation function we fix 
an ion at the origin to find: 

5e/(0,r) = , (2.13) 

Pe 

where p e (r\I) is the (interacting) valence electron density induced by one ion at the origin. A similar relationship 
holds for the ion-ion pair-correlation function, but for the electron-electron pair correlation function the Percus trick 
cannot be used in this form; one cannot "fix" an electron at the orkan. 

The second step follows the basic ideas of the Kohn-Sham schemed, namely that there exists a local single-particle 
external potential w cff (r) which will induce in a non-interacting system the same one-particle density p(r) as is found 
in the full interacting system. This idea holds both for quantum as well as for classical systems. The external effective 
potential felt by species a when species /3 is fixed at the origin follows directly from the Euler equations: 

vf f3 (r)=v a 0(r) + -^j-r-(i™, (2.14) 

where v a p (r) is the direct interaction between species and /i^ x is the excess chemical potential. Thus the electron-ion 
radial distribution function follows from the indirect Kohn-Sham solution of the Euler-equation combined with the 
Percus identity: 

9ei(r) = ! , 2.15) 

Pe 

where Pe(r\v®f) is the density of the unbound (or free) valence electrons obtained from the wave functions that emerge 
from a Kohn-Sham solution o f the .one-center radially symmetric non- interacting Schrodinger equation in the external 
effective potential given by (EyjEi The ion-ion radial distribution function follows from a direct solution of the 



ion-ion Euler equation combined with the Percus-identity: 



9ii{r) = — = ! = exp [-pv II (r)\ , 2.16) 

Pi Pi 



where in the classical context vfJ is commonly referred to as the potential of mean force. Next we expand the 
effective potentia ls (r) and (r) in a functional Taylor expansion around the equilibrium homogeneous densities 
and rewrite Eq. ( 2.14| ) as: 



■Ci(r) = v aI (r) -- p Y.Pij " r'\)h^(r)dr' + ±B aI (r), (2.17) 



where the C Q7 (r) are the homogeneous limits of Eq. (2^), and the Percus trick was used to rewrite (p 7 (r|v 7 /) — p 7 ) 
in terms of the correlation functions h 7 j(r) = <7 7 /(r) — 1. The remaining third and higher order functional derivative 
terms are lumped into the so-called bridge functions B a p{r). We note that these equations do not hold for the effective 
electron-electron potential v ee (r) . 

Up to this point, our formulation is in principle exact. However, since the exact free-energy functionals and the 
related effective external potentials are unknown, some approximations must be made. In the language of the theory 
of classical liquidsEHI, we need a closure relation. For this we follow the approach developed by Chihara, which he 
named the Quantum Hypernetted Chain Approximation (QHNC)E3. The main approximations made by Chihara are 
(roughly in ascending order of importance): 

1. The bare ion-ion potential is taken to be purely Coulombic. This neglects core polarisation effects, but these are 
expected to be small in the metals we study. 

2. The ion-ion bridge function Bji(r) is approximated by the one-component bridge function of an appropriate 
reference state. This is commonly called the RHNC or MHNC approximations, and is expected to be quite 
accurate. We use the repulsive part of the one-component effective pair potential solved in the Percus- Yevick 
approximation as a reference system to calculate the bridge function (see the Appendix for details). 



4 



3. The electron-ion bridge function B e j(r) is set to 0. This is commonly called the hypernetted chain (HNC) 
approximation, and is generally also quite accurate, especially as the electron-ion correlations are expected to 
be weaker than the ion- ion correlations. 

4. The local density approximation (LDA) is used for the one-center electron-ion problem. The calculation of 
the electron-io n co rrelation function reduces to calculating the Schrodinger equation in the external potential 
given by Eq. ( 2.17 ). This is similar to a self-consistent field all-electron calculation for a single atom, except 
that the potential includes not only the nuclear Coulomb contribution, but also terms reflecting the effect of 
the surrounding ions. We solve this effective atomic problem in the LDA, which is widely used in electronic 
structure calculations. The core electrons are treated explicitly, i.e. this is an all-electron calculation. However, 
the core and valence screening effects are separated in a manner similar to the linear unscreening procedure 
used to derive pseudo-potentialsEj. 

5. The valence electron correlations are treated in the jellium approximation. To calculate the full effective po- 
tentials, we need the electron-electron direct correl ation function C ee (fc), which can be re- written in terms of 
the so-called local field factors as was done in Eq. (|2.12 ) where the non-Coulombic correlation part has been 
subsumed into the local field factor G(k). In the QHNC approach, the local field factor is approximated to be 
that of jellium at the average electron density, i.e. it is independent of ionic correlations: 



Cee(fc) = -0Vee(k)[l - G{f {k~ p e )} 



(2.18) 



Thus the electron-electron direct-correlation function uncouples from the other correlation functions in 
Eq. ( |2~IT1) . This approximation is similar in spirit to the LDA approximation and greatly simplifies part of 
the electronic problem, but it is probably the most serious and uncontrolled part of the QHNC closure. 



The appro ximations for the bridge-functions toge ther with Eqns. ( [2.15 ), (2.16), (2.17), and the closure for C ee (k) 



in Eq. ( [2.18 ) reduce the QOZ relations of Eq. ( 2.11 ) to a closed pair of coupled equations for the radial distribution 



functions: 



p e g e i(r) = p e 



(r\v eI (r) - ~ PI J C eI (\r-v'\)h n (r)dr' -^ Pe J C e 



(\r-r'\)h eI (r)dr' 



(2.19) 



g n {r) = expi^ -(3v u (r) - Pl j C n (\v - v'\)h u {r)dv' - p e J C Ie {\v - v'\)h eI {r)dv' + -B{r) ) (2.20) 

which are solved self-consistently. This is the essence of the QHNC approach: the original many-center problem has 
been reduced to an effective one-center problem by replacing the many-body ion-ion correlations with an effective 
external potential that depends self-consistently on the ion- ion correlations. The main advantages are that (a) no 
pseudo-potential is needed, i.e. it is an all-electron calculation and (b) ion-ion and electron-ion correlations emerge 
naturally and on the same footing. Details of the (rather complex) numerical implementation of the QHNC are 
described in the Appendix. 



III. ION-ION AND ELECTRON-ION CORRELATIONS 



A. Ion-Ion radial distribution function 



Armed with the QHNC approach, we can now tackle the electron-ion and ion-ion correlation functions for a set of 
simple metals from the first four rows of the periodic table. As a test of the approach, we compare in Fig. |l| thq-QHNC 
ion-ion radial distribution function for our set of metals to the experimental X-ray data of the Waseda grouped. The 
QHNC provides a faithful representation of gn{r) for all the metals except Ga. The accuracy of the QHNC for the 
other elements suggests that it can also be trusted for Be, for which no experimental data could be found. 

The case of Ga, however, calls for closer examination. While the exact form of gn(r) is sensitive to details of 
the liquid state theory aspects of the closure, i.e. the form of the bridge function, we Jxied various forms of the 
closure without much improvement. On the other hand, when we used the Ortiz-BalloncLil local field factor instead 
of the Ichimaru-Utsumi forjrr£3, a considerable improvement was obtained, in agreement with earlier studies based on 
effective ion-ion potentialstl. This sensitivity of the QHNC approach to details of the local field factor G(q) suggests 
that approximation (5) of the previous section begins to break down. In addition, the d-electrons were very close 
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to being unbound, which made the QHNC algorithm difficult to converge. This instability may be attributed to the 
implicit separation of the exchange-correlation potential into bound and valence contributions in approximation (4), 
i.e. the neglect of non-linear-core-corrections. The fact that the Ortiz-Ballone G(q) seems to work better for Ga is most 
likely-due to an accidental cancellation of errors. It performs considerably worse than the LDA or Ichimaru-Utsumi 
G(g)E3 for the other metals in our set. 



B. Pseudo-atoms and Electron-Ion correlations 



The electron-ion structure factor, defined by Eq. (1.1), can always be re- written in the following fashion: 



S eI (k) = ^Sjj(fc), (3.1) 



which defines a new object, n(k). By taking the Fourier transform we find, using Eq. (1.2), the electron-ion radial- 
distribution function: 

P e9ei{r) = n{r) + p? / n(r - r')gn(r')dr' , (3.2) 



V 

which is proportional to the probability of finding an electron a distance r away from an ion located at the origin. Thus 
a natural interpretation of n(r) is the density of a "pseudo-atom", which, when superimposed according the ion-ion 
radial-distribution function gn(r) gives the correct value of the valence electron distribution. The pseudo-atom is 
independent of ionic correlations only tp.|feKt order in the electron-ion potential, at higher orders it implicitly includes 
3-body and higher order ionic averages™™. 

In the QHNC approximation, the e lectr on-ion radial-distribution function follows directly from the solution of the 
one-body Schrodinger equation (Eq. (2.15)). In Fig. ^| we show these electron- ion radial-distribution functions g e j(xl 



for our set of simple metals. Where possible, they have been compared to ab initio Kohn-ShamO and OF-AIMDEll 
results. As is the case for the ion-ion radial distribution functions, the QHNC approximation gives similar results 
to other methods for all the elements except Ga, where once again an improved agreement is obtained when the 
Ortiz-Ballone G(q) is used. 

It is instructive to compare the pseudo-atom density, included in Fig. || as n(r)/p e , with the electron- ion radial 
distribution function g e i{r). The pseudo-atom density goes to zero for larger r, as it is essentially localised around a 
given ion, while g e i(r) goes to 1 for large r, reflecting the fact that outside the range of the ion's own pseudo-atom, 
g e i{r) simply probes the average density of the pseudo-atoms around the other ions so that the probability of finding 
a valence electron a distance r away is simply related to the probability of finding an ion there. g e i(r) and n(r)/p e are 
essentially identical for small r, as one might expect, while at larger r the effect of the ion-ion weighted superposition 
of the surrounding pseudo-atoms on g e i(f) is evident. 

Because g e i{r) implicitly includes a spherical average, all angular bonding effects are effectively washed out, although 
an indication of the effect of bonding can still be found by comparing geiir) and a superposition of the free-atom 
electron densities!!^!. 

The relationships between the pseudo-atom, the ionic correlations, and the electron- ion correlations become clearer 
in fc-space where the electron-ion structure factor is simply the product of the pseudo-atom density and the ion-ion 



structure factor, as shown in Eq. (3.1). The ion- ion structure factor is sharply peaked at its first maximum k p while 
the pseudo-atom density goes through zero at /co- If < k p , the product form implies that the first peak of S e i(k) 
is negative, and the electron-ion structure is in the so-called low valence class, while ifrifa) > k p , then the first peak of 
S e i{k) is positive, and the electron- ion structure is in the so-called high valence classa^. In Fig. || we plot both the 
electron-ion structure factors S e i(k) and the pseudo-atom densities n(k) for our set of metals. Li, Be, Na, Mg and K 
are in the low-valence class, while Al and Ga straddle the two classes. Only Ca seems to fall outside this taxonomy. 



IV. USING FREE-ATOM FORM FACTORS V.S. METALLIC-ATOM FORM FACTORS 

Neutron scattering probes the fluctuations of the nuclei, while X-ray scattering probes the fluctuations of all the 
electrons. In 1974, Egelstaff etj-qlE3 first suggested exploiting this difference to extract electron-ion correlations for 
liquid metals. In 1987, Chiharaa re-examined the X-ray scattering problem, demonstrating that a careful analysis of 
elastic and inelastic contributions leads to the following coherent scattering intensitycj: 

I x (k) = \f I (k)+n(k)\ 2 S N (k), (4.1) 



G 



where Sjv(&) is the nucleus-nucleus structure factor which emerges, for example, from neutron scattering, fi(k) is 
the ionic form factor, i.e. the ionic electron density, and n(k) is the pseudo-atom density. We shall call the object 
/m(&) = fi{k) + n(k) the metallic-atom form factor. Traditionally the structure factor from X-ray scattering Sx(k) 
has been extracted from scattering intensity as follows: 

Ix(k) = \f A (k)\ 2 S x (k), (4.2) 

where f A {k) is the free-atom form factor, or the free-atom electron-density. 

The difference between the two structure-factors, Spj(k) and Sx(k), stems from the difference between the two form 
factors, f A (k) and = fi(k) +n(k), and provides a measure of the change in electron density upon bonding. In 

Fig. |^ we plot the full free-atom (solid lines), metallic (dashed lines) and ionic (dotted lines) form factors for our set 
of metals. Also included are the pseudo-atom densities (chain lines). The ionic form factor is essentially the same in 
the metallic and the free-atom environments, so the difference between the metallic and free-atom form factors stems 
from the difference between the pseudo-atom density and the valence-electron density of the free atom. 

Because X-rays scatter off all the electrons, not just the valence electrons, the effects of bonding are most pronounced 
when the ratio of the number of valence electrons Z to the total number of electrons Za is high. Thus, as can be seen 
in Fig. [l|, the effects are largest in Li and Be, where the ratios (Z : Za) are (1 : 3) and (1 : 2) respectively, and the 
effect becomes smaller for the other elements, where the ratios are Na: (1 : 11), Mg: (1 : 6), Al: (1 : 4.3), K: (1 : 19), 
Ca: (1 : 10), and Ga: (1 : 10.3). 

In crystalline systems, X-ray studies of charge-densities only provide information on the bonding density for certain 
fixed scattering peaks. Similarly, in liquids the scattering is strongest at the first peak of the structure factor (at 
wave-number k p ), so to observe a difference between Sx(k) and S'iv(fc) it is important that the free-atom and the 
metallic-atom form factors differ near k p . This is demonstrated for Li in Fig. ^ and for Be in Fig. |[ Even though 
the difference between the free-atom and metallic-atom form factors is largest at fc's less than k p , the experimentally 
accessible difference, Sx{k) — SV(fc), is largest at k p . 

In Fig. f^the difference, Sx(k) — Sn(Jc), is shown for our whole set of metals. Generally the peak-height for Sx{k) 
is slightly lower than the peak- height for S]y(k) and the two structure factors are virtually identical away from k p . 
As was anticipated in Ref. ||, the largest difference is for Be, where Sx(k p ) is about 5% lower than Spf(k p ). However, 
Be is extremely toxic, and for that reason its static structure has not yet been measured. Perhaps the best chance of 
observing a difference between Sx(k) and Sxik) is for Li, Mg, or Al, where the difference at k p is about 2%. Another 
possibility includes liquid metallic Si, where the ratio is (1:3.5), and fco is expected to be greater than k p (i.e. Si's 
electron-ion structure is expected to be in the high valence class), so that Sx(k p ) is expected to be larger than SN(k p ) 
and the structure factor may peak in a region where the two form factors differ by a larger amount than is the case 
for the low-valence class metals. 

Measuring these differences will be extremely challenging, since they require two completely different scattering 
techniques, which implies subtracting two different sets of systematic correctioxLSj—.In particular, the removal of 
incoherent scattering effects from the total scattering remains under discussionOoE3. We note that a series of 
experiments measuring the differences between X-ray and neutxon-scattering determinations of Sn(k) have been 
reported for Lo, Na, Mg, Al, Zn, Ga, Sn, Te, Tl, Pb, and Ba. These measurements typically show differences 
that are at best 5 to 10 times laxger than expected from theoretical treatments of the bonding effects, such as those 
shown for the QHNC in Fig. \7j3. In fact for some of the heavier elements, where the Sx(k) — S'jv(fc) is expected 
to be very small due to the large number of core electrons, the differences are severaL-orders of magnitude larger. 
In Fig. 0, we include explicitly the combined X-ray and neutron data of Olbrich et. al£B for Li. Even though their 
differences are smaller than any of the differences measured in the other references cited in [ ^] (in fact they are 
the only measurements which fit within the scale of our graphs£2l), Olbrich et. al. claim that experimental errors 
are too large to see bonding effects in Sn(k). Far .theser^easons, the interpretation of these measurements has 
been called into question by a number of authorsBHaHHLa. The theoretical results are very robust, with simple 
linear response theories in some cases agreeing quantitatively with the much more sophisticated ab-initio Kohn-Sham 
calculations! In a crystalline environment, the Kohn-Sham approach has been shown to agree quantitatively to 
several significant figures with highly accurate experimental measurements of the bonding densities^, suggesting 
that the electron densities calculated within the Kohn-Sham approach for the liquid state analogon of these solid 
state measurements should be highly accurate as well. In fact, for the Kohn-Sham type simulations, finite size and 
statistical finite simulation time effects on the ion-ion structure probably cause larger errors than errors arising from 
the determination of the electron densities. However, these simulation errors are well understood, and will at most 
contribute a few relative % to the difference Sx(k) — SN(k). The considerations above, coupled with the difficulties 
in dealing with the subtraction of two very different sets of systematic corrections to the datacj, lead us to conclude 
that the experiments cited have not yet attained an accuracy sufficient to measure the effects of bonding in liquid 
metals. 
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However, the advent of new high-accuracy X-ray and neutron beam sources coming on line, together with the 
improvement of other techniques such as anomalous X-ray scatteringO, may bring the measurement of these differences 
within experimental reach, at least for a few of the metals in our set. It seems increasingly unlikely that this could 
be measured for many other elements where the ratio Z /Za is smaller and the core-electrons wash out any bonding 
effects. 



V. CONCLUDING REMARKS 



We have carried out QHNC calculations for-Lj—Bc-JNTa, Mg, Al, K, Ca, and Ga. The QHNC formalism, first 
introduced and mainly developed by J. ChiharaEjIiSoEa is a closure to the QOZ relations. Ion-ion and electron-ion 
correlations naturally emerge in a unified fashion, and the interpretation of liquid metals in terms of a "pseudo-atom'' 
helps clarify the meaning of the electron-ion radial distribution functions and structure factors. 

The most serious approximation in the QHNC is probably approximation (5) from section II. B, where the electron- 
electron direct-correlation function C ee (k) is approximated by the form for jellium, making it independent of the 
ion-ion and electron-ion correlations. The sensitivity to the local field factor G ee {k) found for Ga may stem from a 
breakdown of approximation (5), but also from the neglect of non-linear-core-corrections implicit in approximation 
(4). Future work will address both these issues. 

The QHNC reduces to a linear-response formalism if the direct-correlation function C e i (r) / (3 is approximated by 
its low-density or long-range form — w e /(r), suggesting that the., accuracy of the QHNC probably benefits from an 
interference effect which reduces the non-linear response termsoEI. For metallic hydrogen, where the lack of core- 
electrons implies no interference effect, C e i(r)/(3 will differ significantly from its low-density limit. The relative 
importance of non-linear response terms also suggests that approximation (5) may be poor for H. In addition, Xu et. 
al.M showed that small changes in C e j(r)//3 can have a large effect when input into DFT theories of the freezing of 
monatomic H. We expect the DFT theories to be relatively less sensitive to changes in C e i(r)//3 when applied to the 
simple metals in our set. 

The differences between X-ray measurements of the ion-ion structure factor Su(k) interpreted with a free-atom or 
with a metallic-atom form factor are the main experimentally relevant quantities we calculate. This difference, which 
reflects the effects of metallic bonding of the valence electrons, is largest for elements with a large ratio of valence 
to core electrons, such as Li, Be, Mg, Al and maybe Si. To date these bonding effects have not been convincingly 
observed, but with new higher precision instruments coming on line, they may soon be experimentally accessible. 
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VII. APPENDIX: PRACTICAL IMPLEMENTATION OF THE QHNC APPROXIMATION 



A. Overview of the implementation 



In the practical implementation, we follow 2 steps to self-consistency. 

Step 1: the ion-ion loop. For a given g e i(r) and C e i(r), an effective one-component ion-ion effective potential 
is calculated and the 1-component RHNC integral equation is solved self-consistency for gu(r). 

Step 2: the electron-ion loop. For a given gn{r) and the old g e i{r) and C e j(r), an effective electron-ion 
potential v^(r ) is c alculated from Eq. ( 2.17 ). The self-consistent Schrodinger equation is then solved to give a new 
g e i(r) via Eq. (2.15), and the procedure is repeated to obtain self-consistency in g e i(r). 

These two steps are then repeated until full self-consistency is obtained between the two loops. 
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B. Details of the the ion-ion loop 



We first rewrite the ion-ion problem as an effective one-component system with the same radial distribution function: 

g JX (r) = exp[«f?(r)] = g(r) = ex P K ff (r)] (7.1) 



where vj^(r) is the effective potential of mean force for the ions, given by Eq. ( 2.17 ), and wj ff (r) is the effective 
potential of mean force for the one-component system. The equality of the two radial-distribution functions then 
implies that: 

- (3v H (r) + h n (r) - C u (r) + B H (r) = -/3«i(r) + h(r) - C{r) + B(r), (7.2) 

where v\ (r) is the bare potential of the effective one-component system, C (r) is its direct correlation function, and 
Bn(r) and B{r) are the bridge functions of thp_two-component and effective one-component systems respectively. 
We follow Chihara and make the approximations: 

B H {r) » B(r), (7.3) 



which, together with the QOZ relations of Eq. (2.11), implies that: 



/ \ / \ X { ee } (k)\C eI (k)/P\i 
vi (r) = vii {r) M • (7.4) 

Note that if the el ectro n-ion direct correlation function is replaced by its low density or long-range limit: C e i(r)//3 — 
—v e id.r), and Eq. ( 2.18 ) is used for C ee (fc), the effective one-component potential reduces to the usual linear screening 
forma. 

Having now reduced the problem to an effective one-ccpapoaent form (by assuming a fixed g e /(r) and C e /(r)), we 



m 



solve the self-consistent RHNC equations in the usual wayHd'caO. The bridge function is obtained using as a reference 
system the repulsive partpof the Vi(r) solved in the Percus-Yevick approximatiorO which is known to perform well 
for short-range potentials^. The bridge function obtained in this way gives very similar results to the standard 
RHNC approximation (where the reference system is the hard-sphere fluid) for the systems here studied but with the 
advantage that no optimization of the hard-sphere diameter is required. This feature is specially recommended in 
the context of the QHNC theory, where the RHNC equation is solved in combination with the ion-electron integral 
equation. 



C. Details of the electron-ion loop 

Fo r a g iven gn(r) — g(r), and an old g e i{r) and C e /(r), the new effective electron-ion potential follows from 



Eq. ( |2.17D : 



vtf(r)=v e i(r)-p e J Cee(l ^ r ' l) h e i(r)dr> Pl J ^ ^ hn{r)dv\ (7.5) 

where the first term, v e i(r), is the bare electron- ion interaction, the second term describes the screening by the valence 
electrons, both those of the central ion as well as those originating from the pseudo-atom densities of the surrounding 
ions, and the third term describes the interaction of the electrons with the other ions. For the local-field factors 
implicit in C ee (k), we used the Ichimaru-Utsumi forn£3, but exceat for Ga, the simpler LDA form also performed 
quite well. For the bare electron-ion interaction we follow Chiharallj and write: 

vei(r) = + I v ee {\v - v'\)p b e (r)dv> + (i X c[p b e (r) + p B ] - »xc\p e ], (7.6) 



where Z A is the nuclear charge, Pxc[p{ r )\ is the exchange-correlation part of the free energy functional (we take 
the usual LDA parameterisation of Perdew and Zimgcitij of the Ceperley- Alder quantum monte-carlo siraulationa£3), 
and p b (v) is the bound electron density obtained from the solution of the Schrodinger equation. This form is not 
exact within- the LDA, as its derivation implies a linear unscreening process, neglecting the so-called non-linear-core- 
correctionscJ. In fact, this linear unscreening process is not necessary, and the full screening from the combined 
valence and core electron densities can be taken into account, but this will be addressed in a later publication. 







Using this effective electron-ion potential, the one-electron Schrodinger equation: 



2m„ 



t#(r) 



#(r) = c*^(r), 



is solved for the effective potential of Eq. (7.5). The bound electron density is then calculated by means of 



p b e (r\v^)=p b e (r)=J2W b) (r) 

i(b) 



(7.7) 



(7.; 



where the index (b) refers to bound states, while p£(r), the unbound density directly related to g e i{r) through Eq.(2.1!i) 
corresponds to the continuum part of the eigenvalue-spectrum (positive energies) and is calculated as a superposition 
of scattering states. In atomic units this is given byLj& 



pt{r\vf I )=p e + -z / dkk 2 J2W + V[Rli(r)-J?(rk)} 



(7.9) 



where kp is the Fermi wave vector and Rki(r), the radial part of the wave function, is a solution of the equation 

d 2 (rR k i) , [ ;2 1(1 + 1) 



dr 2 



rR M (r)=0. 



Rki{r) must be normalized by its asymptotic limit, i.e. 

lim [rRki(r)] = jf (r k) cos f]i(k) + nf^rk) sin T)i(k), 



(7.10) 



(7.11) 



where ji(x) and ni(x) are_spherical Bessel functions and rji(k) is the phase shift. The phase shifts at the Fermi level 
fulfil the Friedel sum rulgEH: 



-Y / (2l + l)r ]l (k F ) = ZS II (0) 

7T ^ ' 



(7.12) 



where Z is the ionic charge and 5*7/(0) the long-wavelength limit of the ion-ion structure factor. 
We solve the Schrodinger equation in two stages: 

(1) We first lapk for bound states and the eigenvalues of Eq. ( f7~7| ) using the predictor-corrector method on a 
logarithmic griclEil 

(2 ) On ce the bound density is computed (needed to obtain the electron-ion bare interaction via Eq. (|7.6|)), we solve 
Eq. (7.10) for scattering states on a Hcrmann-Skillmann mesh using the Numerov method. Both the logarithmic and 
the Hcrmann-Skillmann grids, as well as the linear mesh in which the correlation functions are stored, span the same 
range of 40.96 a.u. (see below for numerical details). We follow Chihara and introduce a cutoff radius r cut outside 
which the solutions are taken to be of the usual Friedel oscillatory form: 



Vi(r) oc 



47rr 3 



cos(2kpr) 



r > r c 



(7.13) 



where r cut is typically equal to around 4—5 times the electron Wigner-Seitz radius. This procedure avoids difficulties 
with the long-range nature of the potentials while keeping the number of mesh-points needed at a manageable level. 
Typically, we needed 4 or 5 iterations of the ion-ion and electron- ion loops for a maximum tolerance of 0.1% error 
between successive solutions to converge. The screening relation 



J p f e (r)dr = Z 



(7.14) 



is fulfilled with an error around 0.1% for the lower valence metals, and less than 1 %. for ,Ga or Al. Although this is 
somewhat unsatisfactory, and is a larger error than that reported in previous studiesExEj, we found that it does not 
affect the shape of the correlation functions or the effective pair potential at short and intermediate distances. Future 
work will address this issue in more detail. 

We use a linear grid of 4096 points for the ion-ion and electron-ion correlation functions (i.e. a grid size of 0.01 a.u. 
for a maximum distance of 40.96 a.u.). The logarithmic grid contains 983 points whereas the Hermann-Skillmann grid 
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comprises 5 blocks of typically 80, 40, 40, 90, and 800 points respectively. This corresponds to a grid size of 0.0015 
for the first block. The Schrodinger equation for scattering states is solved for 2 5 equally-spaced fc-points between 
and kp and up to 11 plane waves (l m ax — 10). Then, the integral in Eq. (7.9) is computed using Simpson's rule 
with the same number of fc-points. An interpolation by cubic splines is employed to swap between the different grids 
involved in the calculation (linear for correlation functions and logarithmic and Hermann- Skillmann grids to solve 
the Schrodinger equation). Fast Fourier transforms are-used to convert correlation functions from real to reciprocal 
space. In addition, we have implemented Ng's methocE-3 in the ion-ion and the ion-electron loops to accelerate the 
convergence. 



D. Details of the initial setup 

As mentioned before, the implementation of this iterative procedure requires an initial effective pair potential. 
Following Chiharaea, we make use of the so-called jellium-vacancy model (JVM) to obtain such initial potential. The 
JVM can be derived directly from the QHNC approach by the following 2 approximations: 

1. The ion-ion correlation function approximated as a step function: 



h n (r) 



-1 for r < R 
for r>R 



(7.15) 



where R is the ion Wigner-Seitz radius, i.e., R = [3/(47rp/)] 1 / 3 . 



2. The electron-ion DCF used in the effective electron-ion interaction of Eq. (7.5) is approximated to be of a purely 
Coulombic form: 



C e i(r)/P 



Z ie 2 



(7.16) 



These two approximations result in an electron-ion potential v^f(r) which is independent on the ion-ion correlationsEl. 
The Schrodinger equation that follows from this effective potential can then be solved self-consistently as d escribed 
above, and the ensuing g e i( r ) and new C e /(r) used to derive an effective ion-ion pair-potential from Eq. (7.4), for use 
in the initial ion-ion loop. 
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Fig. 1 
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FIG. 1. Ion-ion radial distribution functions gn(r) cal- 
culated by means of the QHNC method (solid lines) and 
compared to X-ray experiments (circles )E2I. The dotted 
lines in the Ga oanel correspond to QHNC with the Or- 
tiz-Ballone G(g)H 
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Fig.2 
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FIG. 2. Electron-ion radial distribution functions 
as obtained from the QHNG-iiftwproximation (solid 
lines), the Orbital- free methodtJiij i{®pen circles) and 
Car-Parrinello Molecular Dynamics£3(crosses). The 
dashed lines represent the pseudo-atom density n(r)/p e . 
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Fig. 3 




k/A~ 1 

FIG. 3. Electron- ion structure factors S e i(k). The 
symbols have the same meaning as in Fig. ^, the dashed 
lines again represent the pseudo-atom density n(k), but 
now in k-space. For the scale, note that n(k = 0) = Z, 
the number of valence electrons. 
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Fig. 4 
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FIG. 4. Free-atom form factors fA{k) (solid lines), 
metallic-atom form factors /m(&) = fi(k)+n(k) (dashed 
lines), and ionic form factors fi(k) (dotted lines), as pre- 
dicted by the QHNC theory. The chain lines represent 
the pseudo-atom density n(k). 
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Fig. 5 
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FIG. 5. The structure factors Sjv(fc) (solid line) and 
Sx(k) (dashed line), of liquid Li. The dotted line cor- 
responds to the difference Sx(k) — S'jv(fc). The metal- 
lic-atom and the free-atom form factor of Li are also in- 
cluded in the figure. 
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Fig. 7 
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FIG. 7. Differences between using a free-atom and a 
metallic-atom form factor to interpret X-ray scattering 
determinations of the static structure factor Su(k) for 
a number of systems as predicted by the QHNC theory 
(note the change in scale from panel to panel). Alterna- 
tively, this can be viewed as the difference between X-ray 
and neutron diffraction determinations of the ion-ion 
structure factor. Also included are the experimental dif- 
ferences between X-ray and neutron diffraction for the 
ion-ion structure factors of Li from reference [ K5| 
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